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Abstract 



We introduce unconditionally stable finite element approximations for anisotropic Allen- 
Cahn and Cahn-Hilliard equations. These equations frequently feature in phase field mod- 
els that appear in materials science. On introducing the novel fully practical finite element 
approximations we prove their stability and demonstrate their applicability with some nu- 
merical results. 

We dedicate this article to the memory of our colleague and friend Christof Eck (1968- 
201 1) in recognition of his fundamental contributions to phase field models. 
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1 Introduction 



The isotropic Cahn-Hilliard equation 

f = V . ( b{u) Vw), w = -eAu + s-^ ^\u) (1) 

was originally introduced to model spinodal decomposition and coarsening phenomena in bi- 
nary alloys, see [[25||27|. Here u is defined to be the difference of the local concentrations of 



the two components of an alloy and hence u is restricted to lie in the interval [—1,1]. More 
recently, the Cahn-Hilliard equation has been used e.g. as a phase field approximation for sharp 
interface evolutions and to study phase transitions and interface dynamics in multiphase fluids, 
see e.g. [[Tj[T6|[23}|46j and the references therein. We note that with 9 = 1 and b{u) = 1 in 
([!]) in the limit 6: ^ 0, we recover the well known sharp interface motion by MuUins-Sekerka, 
whereas 9 = s and b{u) = 1 — leads to surface diffusion; see below for details. 
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The theory of Cahn and HiUiard is based on the following Ginzburg-Landau free energy 

Jn 

where ^ > is a parameter and a measure for the interfacial thickness and Q C M^, = 2, 3, 
is a given domain. The first term in the free energy penalizes large gradients and the second 
term is the homogeneous free energy. In this paper, we consider the so-called zero temperature 
"deep quench" limit, where a possible choice is 

^{u) := J ^ " ""'^ I""! - ^ ' with := f ^/2^s) = | , (2) 
I oc > 1 , J -I 

see [[7|[2T|. Clearly the obstacle potential is not differentiable at ±1. Hence, whenever we 
write "^'{u) in this paper we mean that the expression holds only iov\u\ < 1, and that in general 
a variational inequality needs to be employed. 

We note that ([!]) can be derived from mass balance considerations as a gradient flow for the 
free energy £{u), with the chemical potential w := ^ being the variational derivative of the 
energy £ with respect to u. We remark that evolutions of ([T]) lead to structures consisting of bulk 
regions in which u takes the values ±1, and separating these regions there will be interfacial 
transition layers across which u changes rapidly from one bulk value to the other. With the help 
of formal asymptotics it can be shown that the width of these layers is approximately 6: tt; see 
e.g. [[22[^^. 

In this paper we want to consider an anisotropic variant of £{u), and hence of ([T]). To this 
end, we introduce the anisotropic density function 7 : ^ ]R>o with 7 G C^(]R"\ {0}) n 
C{W^) which is assumed to be absolutely homogeneous of degree one, i.e. 

j{Xp) = \X\j{p) VpgM^VAgM ^ i{p)'P = l{p) Vj9GM^\{0}, (3) 

where 7' is the gradient of 7. Then the anisotropy function defined as 

Aip) = l\j{p)\' VpGM^ (4) 

is absolutely homogeneous of degree two and gives rise to the following anisotropic Ginzburg- 
Landau free energy 

£^{u):= [ l\-f{\/u)\^ + e-^^{u) dx= [ e A{\/ u) + e'^ ^{u) dx ; (5) 
Jn Jn 

see e.g. [[3TJ[35|[37]]. Note that £^ reduces to £ in the isotropic case, i.e. when 7 satisfies 

7(p) = \p\ V p G . (6) 
In this paper, we will only consider smooth and convex anisotropics, i.e. they satisfy 

Hp) . q < l{q) G \ {0} , g G M^ (7) 
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which, on recaUing ([3]), is equivalent to 

7b) + i{p) .{q-p)< l{q) G \ {0} , g G . (8) 

Together with initial and natural boundary conditions, the anisotropic Ginzburg-Landau 
energy ([5]) yields the following anisotropic Cahn-Hilliard equation: 

3v 

e— = V .{h{u)Vw) in Qt := X (0,r) , (9a) 

\c^a-^w = -eV .A'{Vu)+e-^m'{u) in Qt , (9b) 

du Ow 

— = , b(u) — = on dQx (0, T) , (9c) 

u{-,0) = uo in (9d) 

where 9,a ^ ]R>o with a being a factor relating to surface tension in the sharp interface limit, 
and where u is the outer normal to dQ. Moreover, Tio ^ ^ ^ K is some initial data satisfying 
\uo\ < 1. 



An alternative to the no-flux boundary conditions pc^ are the conditions 

Ou 

— = 0, w = g onaf^x(0,T), (10) 
ou 

which are relevant in the modelling of crystal growth. Here in general g G H^{dVt), but for 
simplicity we assume that g = wqq^ G M throughout this paper. 

The anisotropic AUen-Cahn equation based on ([5]) is given by 

e— = e\/ .A'{\/u)-e-^^'{u) in , (Ha) 

yj L 

Ou 

— = 0, on dnx(0,T), (lib) 

u{-,0) = uq in Q . (11c) 



It was shown in [[2}|37| that as 6: ^ the zero level sets of u converge to a sharp interface F 
which moves by anisotropic mean curvature flow, i.e. 

^ V = /.^, (12) 



7(n) 



where V is the velocity of F in the direction of its normal n, and where is the anisotropic 
mean curvature of F with respect to the anisotropic surface energy 



7(n)ds. (13) 

In particular, is defined as the first variation of the above energy, which can be computed as 

:= -V^ . y (n) , 
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i.e. ^ 7(n) = — J^^^^ k.^ V ds; where . is the tangential divergence on F; see e.g. 



Similarly, with the help of formal asymptotics, see e.g. [[T7}|24[|49}[5^ it can be shown that 
the sharp interface limit of (9a-d) with 9 = 1 and b{u) = bo, with 60 ^ ^>o, is given by the 
following MuUins-Sekerka problem 







bo 



dw 



dii 



w 
dw 



Aw 


in fi± , 


-2V 


on r{t) 


a 


on r{t) 





on dQ , 



(14a) 
(14b) 
(14c) 
(14d) 



where Q± denote the domains occupied by the two phases, F = (dQ^) D Q is the interface 
and [-jt denotes the jump across F and with n pointing into the set Of course, if the 
natural boundary conditions ([9c]) are changed to ([10]), then the limiting motion becomes ( |14af - 
c) together with 

w on . (15) 

The problem ( |14a| -c), ( [T5] ) models the supercooling of a molten pure substance, with w playing 
the role of a (rescaled) temperature. Then ( |14b| ) is the so-called Stefan condition, ( |14c| ) is the 
anisotropic Gibbs-Thomson law without kinetic undercooling and ([15]) prescribes the super- 
cooling at the boundary. Here we note that in the quasi- static regime the heat diffusion was 
reduced to Laplace's equation in Q± in ( |14a| ). 

On replacing heat diffusion with particle diffusion, the model ( |14a| -c), ( [T5] ), with w now 
representing a (rescaled) particle concentration, is relevant in isothermal crystal growth where 
a density change occurs at the interfacer, see e.g. []TT|[T3|[T4J and the references therein. We 
remark that in order to recover the anisotropic Gibbs-Thomson law with kinetic undercooling 
in place of ( |14c| ) a viscous Cahn-Hilliard equation needs to be considered, see e.g. [[3|. We will 
look at this in more detail in the forthcoming article p3] |. 

For later use we remark that a solution to ( |14a| -c), ( [T5] ) satisfies the energy identity 

^(^2a J {n) ds- 2 Won yol{n^)^ + bo J |V^|^dx = 0, (16) 

see e.g. fTTl, which is the sharp interface analogue of the corresponding formal phase field 
energy bound 

^(2a — S^{u)-Wdn [ udx]+bo [ \\/ w\^ dx < (17) 

for the anisotropic Cahn-Hilliard equation (9a,b), ( 10) with 9 = 1 and b{u) = bo. 

1 — is given by 



Lastly, the formal asymptotic limit of (9a-d) with 9 = e and b{u) 
anisotropic surface diffusion, i.e. 



V 



(18) 
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where is the Laplace-Beltrami operator on F. The hmit ( fTS] ) in the isotropic case ([6]) was 



formally derived in p6] |, and together with the techniques in e.g. [|37}[38j|34| the anisotropic 
limit ( fTS] ) is easily established. More details on the interpretation of anisotropic sharp interface 
motions as gradient flows for ( [T3] ) and on their phase field equivalents can be found in [|52j]. 

It is the aim of this paper to introduce unconditionally stable finite element approximations 
for the phase field models ([9al-d) and ( |lla| -c). Based on earlier work by the authors in the 
context of the parametric approximation of anisotropic geometric evolution equations [[8}|9|, 
the crucial idea here is to restrict the class of anisotropics under consideration. The special 
structure of the chosen anisotropics can then be exploited to develop discretizations that are 
stable without the need for a regularization parameter and without a restriction on the time step 
size. In particular, the class of anisotropics that we will consider in this paper is given by 



where G 



(19) 



ixd 



for ^ = 1 ^ L, are symmetric and positive definite matrices. We note that 



( [T9] ) corresponds to the special choice r = 1 for the class of anisotropics 



7b) = 



M p G 



r G [1, oc) 



(20) 



which has been considered by the authors in [[9pT|. We remark that anisotropics of the form ( [20| ) 
are always strictly convex norms. In particular, they satisfy ([8]); see Lemma[2[T]below. However, 



despite this seemingly restrictive choice, it is possible with ( [20| ) to model and approximate a 
wide variety of anisotropics that are relevant in materials science. For the sake of brevity, we 
refer to the exemplary Wulff shapes in the authors' previous papers [|8]-[T3|. As we restrict 
ourselves to the class of anisotropics ([19]) in this paper, all of the numerical schemes introduced 
in Section [3} below, will feature only linear equations and linear variational inequalities. The 
numerical approximation of anisotropic phase field models for the class of anisotropics ( [20| ) is 
more involved, and we will consider this in the forthcoming article [15] . 

Let us shortly review previous work on the numerical analysis of discretizations of anisotropic 
AUen-Cahn and Cahn-Hilliard phase field models. Fully explicit and nonlinear semi-implicit 
approximations of the AUen-Cahn equation ( |lla| -c) are discussed in [[31] §8]. In [ [4T| several 
time discretizations for ( |lla| -c) are considered, and unconditional stability is shown for highly 
nonlinear, implicit discretizations. Semi-implicit linearized discretizations are conditionally 
stable on choosing a regularization parameter sufficiently large. Moreover, numerical results 
for anisotropic AUen-Cahn equations have been obtained in e.g. [ |18[|19[|39[|50| |. With particular 
reference to dendritic and crystal growth we mention e.g. p0[[36[|44[[43t [47l|, where a forced 
anisotropic AUen-Cahn equation is coupled to a heat equation for the temperature. We also 
mention the contributions of Christof Eck p2| - [34| |, who introduced homogenization methods 
into the field of crystal growth. 

As far as we are aware, the presented paper includes the first numerical analysis for an 



approximation of the anisotropic Cahn-Hilliard equation (|9a|b). Finally, we mention that nu- 
merical computations for a generalized, sixth order Cahn-Hilliard equation, which is based on 
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a higher order regularization of the energy ([5]) in the case of a non-convex anisotropy density 



function 7, can be found in e.g. [48 , 55 1. 



The remainder of the paper is organized as follows. In Section [2] we consider a stable 
linearization of the gradient A! for the anisotropy function (|4]) and ( [T9] ). This will lay the foun- 
dations for the stable finite element approximations introduced in Section [3| Finally we present 
some numerical results in Section IH 



2 Stable Linearization of A! 



The analysis in this paper is based on the special form ( [T9| ) of 7. Note that for 7 satisfying ( [T9| ) 
it holds that 



N{v)-l{v)l\v). where ^(p) = 5^[7^(p)]-^ G.p VpG]R^\{0}. (21) 



For later use we recall the elementary identity 

2 r (r — 5) = — + (r — sf 



(22) 



Lemma. 2.1. Let 7 he of the form ( |79| ). Then 7 w convex and the anisotropic operator A 
satisfies 



.{p-q)> 7(p) bip) - 7{q)] 

A{p)<'^7{q)J2M<i)]-'h 



VpGM'^\{0},(7G]R^ (23) 
2 Vp G ]R^(7 G E'^V {0}. (24) 



Proof. We first prove (|7]). It follows from ( |2T] ) and a Cauchy-Schwarz inequality that 



7 b) . q = J2^7i{p)]-' {G,p) . q < Y.li{q) = l{q) Vp G \ {0} , g 



1=1 



Together with ^ this implies ([s]), i.e. 7 is convex. Multiplying ^ with 7(p) yields the desired 
result ([23]). Moreover, we have from a Cauchy-Schwarz inequality that 



VpGE^gG]R'^\{0}. 



^=1 



This immediately yields the desired result (|24]), on recalling (|4]). □ 



Our aim now is to replace the highly nonlinear operator A'{p) : — )• in (21 ) with a 



linearized approximation that still maintains the crucial monotonicity property (23 ). It turns out 
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that the natural hnearization is already given in ( [21] ). In particular, we let 

L 

£=1 

L 

L^Ge q = 0. 



B{q) := <^ 



(25) 



Clearly it holds that 



B{p)p = A'{p) VpeM.'^\{0} 



and it turns out that approximating A'{p) with B{q) p maintains the monotonicity property (23 ). 
Lemma. 2.2. Let 7 be of the form ([TP]). Then it holds that 

[B{q)p].{p-q)>j{p) [j{p)-j{q)] Vp,gG]R^. (26) 

Proof. Let p eW^.lf q ^ Oii holds, on recalling ([24]), that 

L L 

[B{q)p] .{p-q)= 7(g) J^^^MT' {p-q)-Gep> 7(g) 5^7^^) {MqT' 7i{p) " 1) 

i=i £=1 

L 

= 7(9) ^^biiQ)]'^ b^{p)? - 7(9) 7b) > 7b) b{p) - 7(9)] . 



If g = 0, on the other hand, then it follows from a Cauchy-Schwarz inequality that 

[B{q)p].{p-q) = [B{q)p].p = Lj2P'GiP = Lj2Mp)f> [j2^i{p)] = b{p)f 

i=i i=i \i=i J 

□ 

Corollary. 2.3. Let 7 be of the form ([TP]). Then it holds that 

[B{q) p].{p-q)> A{p) - A{q) V p , g G . 



(27) 



Proof. The desired result follows immediately from Lemma 2.2 on noting the elementary 
identity ([22]). □ 



3 Finite Element Approximations 

Let {T^}h>o be a family of partitionings of Q into disjoint open simplices a with ha- := 
diam(a) and h := max^-gr'' ^a, so that Q = yj^^-YhO. Associated with is the finite ele- 
ment space 

■= {x e cm : X la is linear V a G T'^} C H\n). 
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We introduce also 

:= {x^S^ \x\ < 1 in f]} C /C := {ry G H^{Q) : |ry| < 1 a.e. in ^} . 

Let J be the set of nodes of and {pj}j^j the coordinates of these nodes. Let {xj}j^J be the 
standard basis functions for that is Xj ^ Xj{Pi) = ^ij ^ We introduce 

TT^ : ^ S'^, the interpolation operator, such that {7r^r]){pj) = r]{pj) for all j G J. A 

discrete semi-inner product on C{Q) is then defined by 



in 



^2)^ := / 7T^{Vi{x)V2{x)) dx 



with the induced discrete semi-norm given by \rj\h := [ {rj, rj)^ ] 2 , for 77 G C{Q). Similarly, we 
denote the L^-inner product over by (•, •) with the corresponding norm given by | • |o. 

In addition to T^, let = to < < • • • < ^a^-i < ^at = T be a partitioning of [0, T] into 
possibly variable time steps := tn — tn-i, n = 1 ^ N. We set r := max^^i^^r r^. 

We then consider the following fully practical, semi-implicit finite element approximation 
for (g-d). Forn > 1 find ([/^, e x such that 
h 

+ (7r^[6([/"-^)] V ly", V x) = V X e 5^ 



I ,X 



(28a) 
(28b) 



where G is an approximation of uq G /C, e.g. = tt^i^q for uq G C(f^). 
Let 

= I s |7(V [/)|^ + (^([/), 1)^ V [/ G (29) 
be the natural discrete analogue of ^ and set 6min := niin5^[_i 1] 6(5). 

Theorem. 3.1. There exists a solution (U'^.W) e x to ^^8^b) with (t/^, 1) = 
([/^~^, 1) ([/°, 1), [/^ Z^" unique. Moreover, it holds that 

^^([/") + (2 a)-^ (7r^[6([/"-^)] V V T^") < . (30) 

/fz addition, ifh^\^ > a^J //"Kt/^, 1)| < J^l dx then is also unique. 



Proof. The existence and uniqueness results follow straightforwardly with the techniques 
in Q, see also [[2TI, on noting from (25) that B{q) G W^^^ is sym metric and positive definite 
for all q G M^. Choosing x = Vl^'' in j28a| ) and x = t^"""^ in ( |28b| yields that 

9 - [/"-\ W)^ + (7r^[6([/"-^)] V VW) = 0, (31a) 

5 (5(V [/^-^) V [/^, V [[/^-^ - [/^]) > (I a-^ W + e-^ [/^-\ [/^"^ - t/^)^ . (31b) 

It follows from ([3Ta[b), on recalling ([22]) and ([27]), that 

\ e |7(V 1^ - \ \UX + Tn{2e a)-' (7r^[6([/"-^)] V V 



lo 2 ^ 



-1 |7-m-l|2 
1^ l/i 



This yields the desired result (30) on adding the constant ^ 1 dx on both sides. □ 
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Remark. 3.2. On replacing (|2Sa]) with 



X 



(32) 



we obtain a finite element approximation for l \lla\-c ). Similarly to Theorem \3. 1 \ existence of a 
unique solution ([/^, W^) G x to (32), {28b), which is unconditionally stable, can then 



be shown. In particular, the solution ([/^, W^) to (32), (28b) satisfies the bound (30) with the 



second term on the left hand side of \30^ replaced by Tn£ ^{\cy\,a ^Y\W^W. 



Remark. 3.3. On replacing the term U^~^ on the right hand side of (28b) with U^, we 
obtain an implicit scheme for which the existence of a unique solution can only be shown 
if the time step satisfies a very severe constraint of the form < C e^9a~^, where the 
constant (7 > depends only on the anisotropy 7 and on the mobility b. In the isotropic case 
^ with constant mobility coefficient h{u) = 60 ^ IK>o Ms constraint can be made precise and 
is given by 

Tn < 2c^£^e{abo)~^] (33) 

see e.g. /[27]/. 



In the remainder of this section we consider the numerical approximation of ( |14a| -c), ( [T5] ). 
In particular, we introduce a finite element approximation for ([9al-c), ([10]). To this end, let 



S^:={xeS^:X = ondn} and S^^ := {x e S"" : x = wen on dQ} . (34) 

We then consider the following fully practical, semi-implicit finite element approximation for 
(l9a|-c), (fTol) with = 1 and b{u) = 60 > 0. Forn > 1 find ([/^, W) e x such that 



Jjn _ jjn-1 



X 



(35a) 



(35b) 



Let 



([/) = 2 a — ^^^([/) - won {U, 1) V [/ G 5^ 



(36) 



Then it holds that the solution to (35a,b) satisfies a discrete analogue to ( 17). 



Theorem. 3.4. There exists a unique solution {U"^, W^) e x to (35a b). Moreover, it 
holds that 

j^!;iu'') + T„ bo I V wX < -^^(f/""') ■ (37) 



Proof. The existence and uniqueness proof is similar to the proof of Theorem |3.1[ but we 
detail it here for the readers' convenience. Let : ^ 
operator for the homogeneous Dirichlet problem on i.e. 



Sq denote the discrete solution 



(38) 
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Hence for [/" G we have that ( 35a ) is equivalent to 

-1 r'hiu"-u"-'^- 



Wan - bo Q 



(39) 



It follows from ( |35bl ) and ^ that f/" e is such that 



> {a-'w9n + e-'U^-\x-U'')'' 



(40) 



where a 



a > 0. There exists a unique G K'^ solving (40) since this is the Euler- 



Lagrange variational inequality of the strictly convex minimization problem 

min { I (5(V t/"-^) V ^^ V z'') + (2 r„ « 5o)"' I V G'^^z'' - t/""') |^ - w^Q + f/""\ z'')''} . 



Therefore, on recalling (39), we have existence of a unique solution ([/^, W^) G x S'^ to 
(|35alb). Choosing x = W/"" - in ([35a|) and x = t^"""^ in ([35b| yields that 



([/^ - [/^-\ - Wdnf + Tn bo (V V T^^) = 



Hence the desired result (37) follows from (22) and (27). □ 



Remark. 3.5. It is easy to show that for U 



n-l 



1 and 



-1 ^ - 



(41) 



^/i^ unique solution to ^35a\ b) is given by = 1 and = wqq. However, if the phase field 
parameter e does not satisfy ( [47] ), then = 1 and wqq is no longer the solution to 

i\35a\ b). In fact, in practice it is observed that for e sufficiently large the solution exhibits a 
boundary layer close to dQ where < 1. This artificial boundary layer, which formally can 
be shown to be also admitted by the continuous problem (i9a\b,d), ([70]), is an undesired effect of 
the phase field approximation for the sharp interface problem i\14a\ -c), ([75]). 



4 Numerical Experiments 

In this section we report on numerical experiments for the proposed finite element approxima- 
tions. For the implementation of the approximations we have used the adaptive finite element 
toolbox ALBERTA, see [51 ]. We employ the adaptive mesh strategy introduced in [16] and [[4|, 
respectively, for = 2 and d = 3. This results in a fine mesh of uniform mesh size hf in- 
side the interfacial region |[/^"^| < 1 and a coarse mesh of uniform mesh size he further away 
from it. Here hf = ^ and he = ^ are given by two integer numbers Nf > N^, where we 

assume from now on that Vt = (—77, 77)"^. As a solution method for the resulting system of 
algebraic equations we use the Uzawa-multigrid iteration from [|4|], which is based on the ideas 



10 



in [ |40| . We remark that recently various alternative solution methods have been proposed, see 
e.g. [[5}^g2H43]]. 



For all the computations we take H = \, unless otherwise stated. Throughout this section 
the initial data G C{Q) is chosen with a well developed interface of width e tt, in which 
varies smoothly. Details of such initial data can be found in e.g. [|4[|6j|T6). Unless otherwise 
stated we always set = 16 tt and Nf = 128, = 16. In addition, we employ uniform time 
steps Tn = T,n=l^N, 



For numerical approximations of ( [T2| ) we employ the scheme ([32]), ( |28b| ). In computations 
for ( 18 ) we use the scheme (28a,b) and fix h{u) = 1 — u'^, 9 = s and a = ^. In all other cases, 
i.e. for the sharp interface limits (14a-c) with ( |14d ) or (15), we fix b{u) = bo = 2, 9 = 1 and 
a = 1 unless otherwise stated. 

For the anisotropics in our numerical results we always choose among 



ANI 



a 

7i{p) = Y.[S'\p\'+p',{l-S')]'^ , with 5>0, 



ANI2: 7 as on the bottom of Figure 3 in Q, 



ANI3 : 7 as on the right of Figure 2 in yjj , 



ANI4: 7 as in Figure 3 in [ 14|. 



We remark that ANI^"^^ is a regularized /i-norm, so that its Wulff shape for S small is given by a 
smoothed square (in 2d) or a smoothed cube (in 3d) with nearly flat sides/facets. Anisotropics 
with such flat sides or facets are called crystalline. Also the choices ANI^, i = 2 ^ 4, represent 
nearly crystalline anisotropics. Here the Wulff shapes are given by a smoothed cylinder, a 
smoothed hexagon and a smoothed hexagonal prism, respectively. Finally, we denote by ANI^ 
the anisotropy ANI^^'^^^ rotated by | in the xi — X2-plane. 



4.1 Numerical results in 2d 



A numerical experiment for ( [T2| ) with the help of the approximation ([32]), ( |28b| ) for the AUen- 
Cahn equation ( |lla| -c) can be seen in Figure [TJ Here the initial profile is given by a circle with 
radius 0.3. We set r = 10"^ and T = 0.05. As expected, the round interface first becomes 
facetted, before is shrinks to a point and disappears. 



A numerical experiment for ( [T8] ) with the help of the approximation ( |28a[ b), for the Cahn- 
Hilliard equation ([9al-d) can be seen in Figure [2} Here the initial profile is given by two circles 
with radii 0.2 and 0.15. We set r = 10"^ and T = 10"^. We observe that the two connected 
components of the inner phase each take on the form of the hexagonal Wulff shape. 

A repeat of the experiment but now for b{u) = 60 = 2, so that the sharp interface limit is 
given by the MuUins-Sekerka problem ( 14a -d), is shown in Figure [3} Here we set r = 10"^ 
and T = 5 X 10"^. Now, in contrast to the evolution in Figure|2| the smaller region shrinks so 
that eventually there is only one connected component of the inner phase. Of course, the final 
interface is converging to the hexagonal Wulff shape. 
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■ ■ ■ 



Figure 2: (ANI3) A phase field approximation for anisotropic surface diffusion ( [T8| ). Snapshots 
of the solution at times t = 0, 2 x 10"^ 5 x 10"^ 10-^ 10"^ A plot of S!" below. 
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0,001 0,002 0,003 0.004 0,005 



Figure 3: (ANI3) A phase field approximation for the anisotropic MuUins-Sekerka problem 
([T4a|-d). Snapshots of the solution at times t = 0, 10"^, 10"^ 4.6 x 10"^ 5 x 10"^ A plot of 
S!^ below. 



The remaining computations in this subsection are for the scheme p5a[ b). In order to visual- 
ize the possible onset of a boundary layer as explained in Remark [33| we present a computation 
for ( 35a ,b) with the initial data = uq = 1. As we set a = 1, the critical value for wqq in 
(pTl) is — ^ = —^IQtt = —64. In our numerical computations this lower bound appears 
tone sharp. In particular, we observe that [/^ = 1 is a steady state whenever wqq > —64, 

^dQ —64 — 10"^. As an example, we present 



IS — — 
3e sharp. 

but a boundary layer forms already for e.g. 
a run for wqq = —65 in Figure |4} where we can clearly see how the boundary layer develops. 
Once the boundary layer has formed, the inner phase first shrinks and then disappears, leading 
to the steady state solution = —1 and = wqq. Note that this phenomenon is completely 
independent from the choice of anisotropy 7. The discretization parameters for this experiment 
were Nf = Nc = 128 and r = 10"^ with T = 10"^ 

Next we simulate the growth of a small crystal, when the sharp interface evolution is given 
by ( |14a| -c), ( [T5] ). In particular, we fix i7 = 8, wqq = —2 and a = 0.03; and we observe that for 
this choice of parameters the condition (41 ) is satisfied if we choose = 32 tt > y vr. A run 
for ( |35a[ b), when the initial seed has radius 0.1, with the discretization parameters Nf = 4096, 
Nc = 128, T = 10"^ and T = 7.5 is shown in Figure [sj We observe that at first the crystal 
seed grows, taking on the form of the Wulff shape of 7. Then the four sides break and become 
nonconvex, with the four side arms that grow at the corners yielding a shape that is well-known 
in the numerical simulation of dendritic growth. 



4.2 Numerical results in 3d 



A numerical experiment for ( [T2| ) in 3d with the help of the approximation ( |32| ), ( |28b| ) for the 
AUen-Cahn equation ( |lla| -c) can be seen in Figure [6} Here the initial profile is given by a 
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Figure 5: (ANI^ ' wgci = —2, e ^ = 32 7t, Q = (—8,8)^) A phase field approximation 
for the anisotropic MuUins-Sekerka problem ( |14a| -c), ( [T5| ). Snapshots of the solution at times 
t = 0, 1, 5, 7, 7.5. A plot ofTi" below. 
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Figure 6: (ANI2) A phase field approximation for the anisotropic mean curvature flow ( [T2| ). 
Snapshots of the solution at times t = 0, 5 x 10"^ 10"^ 1.5 x 10"^ 2 x 10"^. A plot of 8!!j 
below. 

sphere with radius 0.3. We set r = 10~^ and T = 0.03. It can be seen that the initially round 
sphere assumes the cylindrical Wulff shape as it shrinks, before the interface shrinks to a point 
and disappears completely. 



A numerical experiment for ( [T8] ) with the help of the approximation ( |28a[ b) for the Cahn- 
Hilliard equation ([9^d) can be seen in Figure |7| Here the initial profile is given by a sphere 
with radius 0.3. We set r = 10"^ and T = 10"^. We can clearly see the evolution from the 
round sphere to the strongly facetted Wulff shape. 



A numerical approximation for the sharp interface problem ( |14a| -d), is shown in Figure [8j 
Here the initial interface is given by the boundary of a 8 x 1 x 1 cuboid with minor side length 
0.1, and we set r = 10"^ and T = 5 x 10~^. We can observe that during the evolution the 
elongated facets become bent and nonconvex, before the solution converges to the Wulff shape. 
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Figure 7: (ANI4) A phase field approximation for anisotropic surface diffusion ( fTS] ). Snapshots 
of the solution at times t = 0, 10-^ 2 x 10-^ 5 x 10-^ IQ-^. A plot of below. 
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